Breaking the exponential wall in classical simulations of fidelity 
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We analyze the efficiency of available algorithms for the simulation of classical fidelity and show that their 
computational costs increase exponentially with the number of degrees of freedom for almost all initial states. 
Then we present an algorithm whose cost is independent of the system's dimensionality and show that, within a 
continuous family of algorithms, our algorithm is the only one with this property. Simultaneously we propose a 
general analytical approach to estimate efficiency of trajectory-based methods. 
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Introduction. While the solution of the time-dependent 
Schrodinger equation scales exponentially with dimension- 
ality and is feasible for only a few continuous degrees of 
freedom (DOF), classical (CL) molecular dynamics simula- 
tions are, in principle, feasible for millions of atoms. It may 
therefore be surprising that papers studying classical fidelity 
(CF) have provided numerical results for only one or a few 
DOF fl'-'T]. A notable exception is Ref. [5], which, for the 
largest systems, relies on initial densities given by character- 
istic functions. Below we explain this situation by showing 
that not only quantum (QM) but also all previously used CL 
algorithms for fidelity scale exponentially with the number D 
of DOF for initial states other than characteristic functions. 
Hence even when QM effects are negligible and CL picture is 
appropriate, the "simple" CL simulations may be unfeasible. 
Since numerical simulations are important for testing analyti- 
cal theories of CF in large systems, we design an efficient CF 
algorithm that avoids the exponential scaling with D. 

Quantum and classical fidelity. While important in its own 
right, CF can be viewed as the CL limit of quantum fidelity 
(QF) [6], introduced by Peres [7J to measure the stability of 
QM dynamics (QD). QF is the squared overlap Fqm [t) at time 
t of two quantum states, identical at i = 0, but evolved with 
two different Hamiltonians, Hq and = Hq + eV: 

^QmW I/qmWI', (1) 

/qmW :=(V'|t/r*f^o|V'), (2) 

where /qm(^) is the fidelity amplitude and [/* :— 
exp(— iiJft/fi,) the QM evolution operator. Rewriting Eq. ^ 
as /qm(0 = (V' V') with the echo operator [/* := U~*Uq, 
it can be interpreted as the Loschmidt echo, i.e., an overlap 
of an initial state with a state evolved for time t with Hq 
and subsequently for time —t with H^. (In general, we write 
time i as a superscript. Subscript e denotes that was used 
for dynamics. If an evolution operator, phase space coordi- 
nate, or density lacks a subscript, Loschmidt echo dynamics 
is implied.) QF amplitude (|2j is ubiquitous in applications: it 
appears in NMR spin echo experiments 11811, neutron scatter- 
ing 19], ultrafast electronic spectroscopy OlOll, et c. QF ([T]l is 
relevant in QM computation and decoherence 111 ill , and can be 
used to measure nonadiabaticity 1I12I1 or accuracy of molecular 
QD on an approximate potential energy surface ifisll . 



Definition O can be generalized to mixed states in differ- 
ent ways 10, [l3l, but we will assume that the initial states 
are pure. In this case, one may always write QF ^ as 
FQuit) = Tr (pIpq) where p* := U*pU~* is the density op- 
erator at time t. In the phase-space formulation of QM me- 
chanics, QF becomes i^QM(i) = h^^ j dxp\y^{x)pQ^{x) 
where x {q,p) is a point in phase space and Aw(a;) := 
/ di{q-il2 \A\ g + C/2)e^P«/^is the Wigner transform of A 
This alternative form of QF provides a direct connection to its 
CL Umit, which is precisely the CF, defined as HI 01 



i^CL(i) F^it) = h 

= -Fecho(0 = h 
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where the first and second line express CF in the fidelity and 
Loschmidt echo pictures, respectively, is the CL phase- 
space density evolved with H^, and is this density evolved 
under the echo dynamics. We omit subscript "CL" for CL 
quantities F and p since CF is the main subject of this paper 
Algorithms. The exponential scaling of QD with D is well 
known. As for CF, Eqs. ©-(HI) may be evaluated, e.g., with 
trajectory, grid, or mesh-based methods. Clearly, the grid- 
based methods would suffer from a similar exponential scal- 
ing as QD on a grid. We focus on the most general and 
straightforward trajectory-based methods, which are obtained 
from Eqs. ©-dill using the Liouville theorem, yielding equiv- 
alent expressions 



dx^p{x^ ^)p{Xf^ *) and 



i^echo(i) - h-"" / dx''p{x-')p{x''). 
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:= <I>*(a;°) where is the Hamiltonian flow of 

where := $7* o $^ is the Loschmidt 
echo flow. Since it is the phase space points rather than the 
densities that evolve in expressions (HI)-®, we can take p = 
pw, i.e., the Wigner transform of the initial QM state. We 
further rewrite Eqs. ©-(ISll in a form suitable for Monte Carlo 
evaluation, i.e., as an average 



Above, x* 
and a;* :— 
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where W is the sampling weight for initial conditions a-". 
The weight can be any positive definite function, but it is ad- 
vantageous to consider the weight to be related to the den- 
sityp. While previously used algorithms sampled from p 
IS |4j,|5|], we consider more general weights — Wi^[{x^) := 
p{x°)'^ and W = Wm{xo*) = p($o *(a;°))*^ for the echo 
and fidelity dynamics, respectively. These weights yield M- 
dependent algorithms 



FM^M{t) = iM{p{x;')p{x^'y-^' 
K,i,,.Mit) ^ iM{p{x-')p{xy-^' 



-*)M, (7) 

p(^0)M, (8) 

where Im '■— J p{x^)'^'^ dx^ is a normalization factor In 
both families of algorithms Q-®, sampling can be done by 
Metropolis Monte Carlo for general dynamics and any pos- 
itive definite weight p*^. For M > 0, the echo algorithms 
(|8]l are, however, much more practical since the initial state 
is often known explicitly (and generally is much smoother 
than the final state), making sampling easier. Furthermore, for 
simple initial states such as Gaussian wavepackets (GWPs), 
the Metropolis sampling in the echo algorithms can be re- 
placed by analytical sampling. Therefore, for M > the 
fidelity algorithms are more of a theoretical possibility than 
a practical tool. For M = 0, the sampling is uniform and 
makes sense only for a compact phase space of finite volume 
fi = i7f = [nih)^ where Vli and ni are respectively the 
phase space volume and Hilbert-space dimension for a single 
DOF. For M > 0, importance sampling based on the weight 
Wm is used and an infinite phase space is allowed. For gen- 
eral M, the sampling is only defined for CL states (such as 
GWPs), for which p > 0. However, for M — and for the 
important special case of M — 2, the sampling is defined for 
any pure state, i.e., even for negative values of p. 

In order to compute CF directly from algorithms (|7]i or 
dHJ, the normalization factor Im must be known analytically. 
For general pure states, Im is known analytically only for 
M = 0, 1, or 2. For M = 0, /q = nf because of the 
requirement of finite phase space. For both M = 1 and 
M = 2, hi = 1 since Trp = Tr p^ = 1. For M ^ 
{0, 1, 2}, algorithms (|7]l and (O can only be used for special 
initial states. E.g., for initial GWPs /9(a;) = g{x;X,a) := 
2-°exp[-(q-Q)Va2-(p-P)2aVS2] where X is the 
center and a the width of the GWP, we have Im = 
(2*^^^/M) for general AI > 0. However, the unknown 
normaUzation factor can be removed from Eqs. dTjl and (HJ 
by dividing them by the value of I2 [note that hiO) ~ ^(0)] 
obtained with the same algorithm and trajectories. Resulting 
"normalized" (N) algorithms. 
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are practical for general initial states and for any M. As far as 
we know, from the four families of algorithms dTji, ([8]), (|9]l, and 
( [Tol l only echo-1 ([8]l has been used previously [2, 4, 5]. Note 
however, that for initial states given by characteristic func- 
tions, echo-1 = echo- A/ = echo-N-i\/ for all M > 0. 

Efficiency. The cost of a typical method propagating N tra- 
jectories for time t is 0{c{tN) where Cf is the cost of a sin- 
gle force evaluation. However, among the above mentioned 
algorithms, this is only true for the fidelity algorithms with 
M = 0. Remarkably, in all other cases, the cost is 0{cft^N). 
For a single time t, the cost is linear in time, but if one wants 
to know CF for all times up to t, the cost is quadratic with 
t. For the echo algorithms, it is because one must make full 
backward propagation for each time between and t. For the 
fidelity algorithms, it is because the weight function p{x^*')^'^ 
changes with time and the sampling has to be redone from 
scratch for each time between and t. In other words, differ- 
ent trajectories are used for each time between and t. 

The above estimates are correct but not the full story. There 
are hidden costs since the number of trajectories N required 
for convergence can depend on D, t, dynamics, initial state, 
and method. One usually empirically increases N until con- 
vergence, but this is often impracticable. Instead, we esti- 
mate N analytically. An essential point is that N is fully 
determined by the desired discretization error CTdisci - The ex- 
pected systematic component of (Tdisci is zero or 0{N~^) for 
all cases studied and is negligible to the expected statistical 
component a = 0{N^^/^) which therefore determines con- 
vergence. This statistical error is computed as a'^{t,N) = 
2 

F{t, N)'^ — F{t, N) where the overline denotes an average 
over infinitely many independent simulations with N trajec- 
tories. Hence we can formulate the problem of efficiency pre- 
cisely: "What N is required to converge fidelity F to within 
a statistical error cr?" We let be a function of F because 
in many applications, one is interested in F above a certain 
threshold value F„^i„. This threshold can vary with applica- 
tion: it may be close to unity (in quantum computing) or to 
zero (yet finite, in calculations of spectra), but in general will 
be independent of D. 

The discretized form of Eq. (|7]i is FM-M{t,N) — 



ImN ^J2f= i PcL(2:,j)pcL(a;oj) 
FM-M{t,Nr - ll,N-\p{x:yp{x^y-^''] 
(1 - iV-i)F2. Similarly, from Eq. ^ Fecho-M(t, A^) = 
lMN-^j:f^,pixj')p{x''^)^-'', hence Fecho-M(t, iV)^ = 
/|,iV-i(p(:r-*) V(a:°)^-^^Op (.o) M + (1 - N- ^)F\ 

Realizing that Ffid^ A/ (t,A^) = ^;cho-M(t, A^) F{t) in 
both cases, we obtain the same error 



t\l-M 



from 



which 



Jm 



echo- A/ 



= N-\ImJm-F^), 



dx"" p{x-'f pix") 



0\2-M 



In the special case of M — 2, we find our main result, 



(10) 



^fid-2 — '^echo-2 



(1~F^) 



(11) 

(12) 



(13) 
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This expression shows that for general states and for general 
dynamics, statistical error of Ffid-2 or of Fecho 2 depends only 
on N and F. In other words, the number of trajectories needed 
for convergence is independent of t, D, or dynamics of the 
system. This important result is due to the fact that for the 
sampling weight W = p^, each numerical trajectory con- 
tributes evenly to the weighted average (at time t ~ 0). 

As for algorithms (|7]i-(|8]l with M 7^ 2, one might hope to 
improve convergence by employing the normalized versions 
(l9Tl-(fT0]i. The error analysis is simplified using the formula for 
statistical error of a ratio of two random variables. 



0"B 



AB- AB 



AB 



(14) 



In our case, F^^uit,^) AjB where A = FM{t,N), 
B = Fm{0,N), a = F{t), B = F{0) = 1, and a a 
and aB are given by Eq. (fTTT i. The only unknown in 
Eq. ( fT4] i is AB. For the normalized echo algorithms 
(doll, we have AB = Feeho-M(i, iV)^;cho-A/(0, A^) - 
ll,N-\pix-')p{x'>r-^'^),^,o.M+{l~N-')F{t)F{0)^ 
N-^ImKm + (1 - N'^)F where Km 



J dx^ p{x *)p{x 



0\3-M 



The same derivation goes 



through for the fidelity algorithms. The final error can in both 
cases be written as 



= ( Jm - 2KmF + U^mF^ 



(15) 



Exponential growth of the error for M ^ 2. Now we will 
show that the special case M = 2 is unique and that all the 
other above-mentioned algorithms (which include all the al- 
gorithms available in the literature) have an error growing ex- 
ponentially with D. Since we are searching for counterexam- 
ples, special cases are sufficient. For us these will be initial 
GWP states and "pure displacement" (PD) or "pure squeez- 
ing" (PS) dynamics [3]. All calculations can be done analyti- 
cally using special cases of the integral 

J dqexp[-ci{q - qif - C2{q - ^2)^] 



Cl + C2 



D/2 



exp 



C1C2 , s2 



Cl + C2 



In the PD case, the center of the GWP moves while 
both its shape and size remain constant. Such fidelity 
dynamics can be realized exactly by two displaced sim- 
ple harmonic oscillator (SHO) potentials with equal force 



constants. 



For PD, the width af, = 



= a = 



a 







= a and either X* = + AX* or X* = X° 



AX*. CF is F{t) 



exp 



J dx g{x; X*^ , a)g{x; X'^ ,a) — 
and the factor (fT2] l 



f AP*a 
{ ^ 

needed in the statistical error can be expressed in terms of F 

as Jm = ( ) F''-' with 7m = 4 - 8/ (4 - M) . Using 



this result in Eq. ( ITTT l. the statistical errors are 



2 2 

'^fid-A/, PD ~ '''echo-Af . PD 



N 



/3o = 2ni and /3o<a/<4 



{PZf-"^' - F') 
4 



(4 - M)M 



(16) 
(17) 



Note that Pm > 1 and the minimum /32 = 1 is achieved for 
M = 2. The minimum agrees precisely with the general result 
(flJl l. Except for M — 2, (Sm > 1, showing that even in the 
simple case of PD dynamics, the errors of all algorithms from 
the families (|7]i and dD grow exponentially with D, which is 
the second major result of this paper. The normalized methods 
(|9]l and (fTOl l lower the prefactor of the error but do not change 
the exponential scaling with D: Since Km = [2'^^^^/(4 — 
TV/)]!?^*^/-! where 6m = 3 - 2/(4 - M), statistical errors 



"'fid-N-A/, PD 



''echo-N-Af, PD 



= N^'^^M {F'""' - 2F^ 



In the PS case, the center of the GWP remains fixed while 
its width narrows in some directions and spreads in others. 
Such fidelity dynamics is realized exactly by two inverted 
SHOs with common centers and different force constants. An- 
alytical calculations show that the errors of different algo- 
rithms again grow exponentially with D (see Table I). 

To summarize, in all cases studied, for D ^ 1 the number 
of trajectories required for a specified convergence is 



N = (T-^a{F)P^ 



(18) 



where a and /3 depend on the method and dynamics and are 
listed in Table J] For both fidelity and echo algorithms with 
M — 2, for any dynamics and any initial state, the coefficient 
13 = 1, implying independence of D. Note also that algo- 
rithms with M = 2 are automatically normalized. For all 
other algorithms (both echo and fidelity, both unnormalized 
and normalized, and for any M 7^ 2) and for both PD and PS 
dynamics, /3 > 1, implying an exponential growth with D. 
This growth is dramatic for M = (/3 = 2ni ^ 1): since nf 
is the Hilbert space dimension, the cost of M = algorithms 
approaches that of QF. This is unfortunate since Ffid-o is the 
only algorithm that scales linearly in time. On the other hand, 
for the most intuitive and most common M = 1 algorithms, 
/3 = 4/3 or \/2, and the growth is much slower, although still 
exponential. We cannot exclude existence of a faster CL al- 
gorithm; however, we doubt existence of a CL algorithm that 
would be both linear in t and independent of D. 

Numerical results and conclusion.To illustrate the analyti- 
cal results obtained above, numerical tests were performed in 
multidimensional systems of uncoupled displaced SHOs (for 
PD dynamics), inverted SHOs (for PS dynamics), and per- 
turbed kicked rotators (for nonlinear integrable and chaotic 
dynamics). The last model is defined, mod(27r), by the map 
qj+i = qj+Pj,Pj+i = pj -VW (qj+i)- eVV{qj+i) whem 
W{q) = —kcosq is the potential and V{q) = — cos(2(7) 
the perturbation of the system; k and e determine the type of 
dynamics and perturbation strength, respectively. Uncoupled 
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TABLE I. The number of trajectories needed to achieve a given a is 
for D > 1 given by iV = a-^a{F)l3°. The table Usts a{F) and 
j3 for different cases. Note that fid-0, echo-1, echo-l', and echo-N-1 
results are for initial GWPs and exhibit exponential scaling with D 
while echo-2 result, valid for any initial state, is independent of D. 
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FIG. 1. Convergence of different classical fidelity algorithms 
in a 100-dimensional system of perturbed (e — 10^*) quasi- 
integrable (k — 0.2) kicked rotators with ni — 131072. Algo- 
rithm echo-2 agrees with the QM result and converges with only 
— 2048 trajectories whereas the echo-1, echo-l', and echo-N- 
1 results are far from converged even with A ~ 7 x 10^. Fully 
converged CL{N = oo) is computed as a product of 100 one- 
dimensional fidelities. The "hopelessly" unconverged fid-0 algo- 
rithm not shown. For clarity, echo- 1 ' error bars not shown for t > 20. 



systems were used in order to make QF calculations feasible 
(as a product of D 1 -dimensional calculations); however, the 
CF calculations were performed as for a truly ZJ-dimensional 
system. The initial state was always a multidimensional GWP. 
Expected statistical errors were estimated by averaging actual 
statistical errors over 100 different sets of N trajectories. No 
fitting was used in any of the figures, yet all numerical results 
agree with the analytical estimates. Note that Table I and fig- 
ures show results for algorithm echo-l', 

^;cho-r(t) = 1 + {p{x~') - P(^°))p(,0), 

which is a variant of echo-1 accurate for high fidelity. Both 
echo-1 and echo-l' reduce to echo-N-1 if normalized. 

Figure [T] displays fidelity in a 100-dimensional system of 
kicked rotators. It shows that echo-2 converges with several 
orders of magnitude fewer trajectories than the echo-1, echo- 
r, and echo-N-1 algorithms. Figures |2] and [3] confirm that 
o-echo-2 is independent of D while (Techo-i, cTecho-i', and decho-N-i 
grow exponentially with D. The normalized echo-N-1 algo- 
rithm is the most efficient among the methods with M = 1. 

To conclude, we have shown that not only QF, but also CF 
algorithms can be unfeasible in complex systems due to the 
exponential scaling with dimensionality. We have proposed 
an efficient CF algorithm for which this exponential scaling 
disappears. In the special case of initial densities given by 
characteristic functions all echo-Af and echo-N-A/ algorithms 
(for M > 0) collapse into a single algorithm. In partic- 
ular, the "natural" algorithm sampling from p is equivalent 
to our algorithm sampling from p^. This may explain why 
high-dimensional calculations were previously done only with 
characteristic functions. These results should be also useful in 
applications computing more general overlaps of phase space 
distributions. Finally, we have described a technique to an- 
alyze efficiency of general trajectory-based algorithms. This 
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FIG. 2. Statistical error grows exponentially with D for the echo- 
1, echo-l', and echo-N-1 algorithms and is independent of D for the 
echo-2 algorithm. Dynamics corresponds to pure displacement, A « 
10^, and time was chosen separately for each D so that F ~ 0.3. 
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FIG. 3. Statistical error grows exponentially with D for the echo-1, 
echo-l', and echo-N-1 algorithms and is independent of D for the 
echo-2 algorithm. Dynamics corresponds to pure squeezing, A ~ 
10*^, and time was chosen separately for each D so that F ~ 0.99. 
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can be useful in developing approximate methods for QD of 
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